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We describe the formation and properties of Holstein polarons in the entire parameter regime. Our 
presentation focuses on the polaron mass and radius, which we obtain with an improved numerical 
technique. It is based on the combination of variational exact diagonalization with an improved 
construction of phonon states, providing results even for the strong coupling adiabatic regime. In 
particular we can describe the formation of large and heavy adiabatic polarons. A comparison of the 
polaron mass for the one and three dimensional situation explains how the different properties in 
the static oscillator limit determine the behavior in the adiabatic regime. The transport properties 
of large and small polarons are characterized by the f-sum rule and the optical conductivity. Our 
calculations are approximation-free and have negligible numerical error. This allows us to give a 
conclusive and impartial description of polaron formation. We finally discuss the implications of our 
results for situations beyond the Holstein model. 

PACS numbers: 71.38.-k,71.38.Cn,71.38.Ht 
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I. INTRODUCTION 

In many materials it is the strong coupling of elec- 
tronic and lattice degrees of freedom that determines the 
charge transport or the optical properties^—. A basic ef- 
fect is the renormalization of the charge carrier mass 7 -. 
The mass renormalization is especially large in materials 
with a short-ranged interaction of electrons with optical 
phonons, as in molecular crystals or the colossal mag- 
nctorcsistivc manganitcs 8 -. In such materials it happens 
that charge carriers, which otherwise resemble free elec- 
trons apart from a certain mass renormalization, give 
way to a kind of new quasiparticle with entirely different 
properties, the Holstein polaroid. The mass of Holstein 
polarons can exceed that of a Bloch electron by several 
orders of magnitude, and may thus explain a low conduc- 
tivity or thermally activated transport^. On the other 
hand, polarons have specific optical signatures^. 

From early work we understand the fundamental mech- 
anism of polaron formation very well in two opposite 
cases: In the antiadiabatic regime of large phonon fre- 
quencies, where strong coupling perturbation theory pro- 
vides the picture of a Lang-Firsov polaroid, and in the 
strict adiabatic limit of static oscillators, where the elec- 
tron and phonon wavefunctions facto r 13 ' 14 . 

Much analytical and numerical work has been per- 
formed to extend our understanding to finite phonon 
frequencies and the intermediate regime, including ex- 
tended variational approaches^—, exact diagonaliza- 
tion (ED)i£~— , density matrix renormalization group2£, 
quantum Monte Carlo^— , and different perturbative or 
approximate schemes 3 ^— . Numerical calculations pro- 
vide reliable data for a wide range of phonon frequencies 
and coupling strengths, against which any discussion of 
Holstein polarons must be checked (for a recent review, 
see e.g. Refs. I36ll37l ). It is however very difficult to per- 
form comparable calculations for small phonon frequen- 



cies, since the increasing number of phonon excitations 
poses a severe restriction on any, even large-scale, nu- 
merics. So far this prevented a thorough study of the 
experimentally most relevant adiabatic regime, in partic- 
ular of large and heavy polarons as predicted by adia- 
batic approximations for one-dimensional (ID) systems. 
Also, numerical data for optical or transport properties 
are scarc o 38 ' 39 . 

In the present contribution we use a new numerical 
technique to extend previous studies deep into the adia- 
batic regime. The key idea is an improved construction 
of phonon states that allows for the simultaneous rep- 
resentation of small and large oscillator shifts with few 
elements in a Hilbert space. The construction overcomes 
the problem that the number of oscillator energy eigen- 
states needed grows fast at small frequency. We com- 
bine this construction with the phonon state selection 
in variational exact diagonalization (VED)^. Based on 
this method we present here results for phonon frequen- 
cies down to a hundredth of the electron transfer integral 
and still at large coupling, extending the frequency range 
accessed in previous studies by more than one order of 
magnitude. We explore the Holstein polaron in parame- 
ter regimes that, so far, were only accessible to approxi- 
mate treatments. The validity of common concepts, and 
thus our very picture of polaron formation, can now be 
assessed. This includes the question which properties of 
the static limit remain at small but finite phonon fre- 
quencies, and how the limit itself is approached. To our 
knowledge it is for the first time that the formation of 
large and heavy ID polarons is studied with unbiased 
numerics. 

The paper is organized as follows. In Sec. [H] we explain 
our numerical technique in application to the Holstein 
model. The differences between the ID and 3D case are 
discussed through the comparison of the polaron mass in 
Sec. IIII1 and the polaron radius in Sec. IIV1 The study of 
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the ID case is continued with the f-sum rule and optical 
conductivity in Sec. [V] We conclude in Sec. IVIII with a 
discussion of our results in the broader context of polaron 
formation under general conditions. 

II. THE HOLSTEIN MODEL AND THE 
NUMERICAL METHOD 

The Hamiltonian of the Holstcin model on a chain (cu- 
bic lattice) in dimension D = 1 (D = 3) is given by 

D 

H = —t c jL+ ei C n + C n C n+ei 

f=l n (1) 

~ \/epUaY^( b l + K) c l c n + w b n b n ■ 
ri n 

The fermionic operators cj, create the electron, and the 
bosonic operators 6J, create a phonon at site n. The 
sums run over all lattice sites n, with the vectors along 
the lattice axes. The Holstein model has three param- 
eters: the electron transfer integral (or hopping matrix 
element) t, the phonon frequency luq, and a third pa- 
rameter e p , A = e p /2Dt, or g 2 = e p /ujQ specifying the 
electron-phonon coupling strength. In the following we 
give all numerical values with the hopping matrix ele- 
ment as the unit of energy, so that t = 1. All results are 
for a single electron at zero temperature. 

Our calculations for the Holstcin model are performed 
with standard ED methods, such as the Lanczos algo- 
rithm and the Kernel Polynomial Method^. The fun- 
damental task is the construction of an appropriate sub- 
space of the complete infinite Hilbert space of the single 
fermionic and many bosonic degrees of freedom of the 
Hamiltonian ([1} . A very efficient construction is provided 
by VED, which achieves extreme accuracy not only for 
the polaron groundstate, especially in the regime of in- 
termediate coupling and phonon frequency^, but also for 
models with different types of fermion-boson coupling^. 

VED is based on an increasing sequence of subspaces 
of the complete Hilbert space. A certain subspace con- 
tains all basis states that are produced through a given 
number A gcn of applications of the Hamiltonian to the 
phonon vacuum, that is either by excitation of a phonon 
or an electron move. For example, it contains the state 
with iVgen phonons at the electron position, and the one 
with a single phonon N gcn — 1 lattice sites away from the 
electron. The number of elements in the subspaces grow 
only by a factor close to D + 1 as -/Vg en increases by 1 , and 
already small subspaces provide a good approximation of 
the true groundstate. All quantities, e.g. variational es- 
timates of the groundstate energy, can be improved by 
increasing N gen until convergence. 

The main restriction of the VED construction, com- 
mon to most ED studies, is the large number of phonons 
at large couplings or small phonon frequencies. To 
overcome this restriction we will replace the standard 
Fock basis with a more general construction. For a 



single harmonic oscillator (everything generalizes imme- 
diately to multiple oscillators at different lattice sites) 
the Fock states with fixed boson number are the states 
\n) = (&t) n /\/?i!|vac). If a force, say — ^/e p w (6^ + b), is 
applied the average number of bosons in the groundstate 
is given by g 2 = e p /ujo. The larger e p or the smaller 
ljq , the more bosons are needed to account for the os- 
cillator elongation (b' + b) = 2g. On the other hand, 
we know that the groundstate is simply a coherent state 
\g) c = S(g)\v&c), obtained from the vacuum with the 
shift operator S(a) = cxp(a&^ — a*b). 

In the polaron problem the effective force exerted on 
the oscillator, hence its state, changes as the electron 
moves. We therefore should try to find states covering the 
entire range of oscillator elongations. This motivates the 
following construction of shifted oscillator states (SOS): 
To a given parameter a > 0, start with the coherent 
states \na) c = S(na)\vac). Orthogonalize these states 
with the Gram-Schmidt procedure, in the given order, to 
obtain the SOS \n) a . The state |0) CT = |0) c remains the 
vacuum state, while the higher states are a complicated 
mixture of coherent states. The Fock states are recovered 
in the limit a -> 0, when S(a) ~ 1 + a(tf — b). More 
properties of the SOS are given in the appendix (Sec.lA")). 

The average elongation and energy of the SOS \n) a 
grows with n, so that they can replace the Fock states in 
the VED. The optimal parameter <j can be determined 
variationally through minimization of the groundstate 
energy. For our purposes, the choice a = g/N gen is close 
to optimal, and was used in all calculations. For larger 
coupling g, hence larger optimal a, the SOS account for 
larger oscillator elongations. Conversely, if the maximal 
number of states iVg en per oscillator is increased and a de- 
creases, the SOS sample ever finer details of the oscillator 
wavefunctions. In fact, since in our calculations always 
cr < 1, the states \na) c overlap in position space, and the 
SOS do not resemble coherent states after orthogonaliza- 
tion. Instead they interpolate between 'quantum' Fock 
states generated by w and 'classical' coherent states gen- 
erated by S(a). Note that a construction with a fixed 
small number of coherent states would have to fail: For 
small ujq it is not possible to represent all oscillator wave- 
functions - with small width and different positions - if 
only few states are available. For the SOS the necessary 
increase in accuracy is accomplished by construction if 

-Ngen grOWS. 

For the polaron problem, the gain in efficiency ob- 
tained using SOS is demonstrated in Tab. E Only this 
efficiency gain allows us to perform calculations for fre- 
quencies ujq <C 0.1 at A = I, extending the accessi- 
ble parameter range by more than one order of magni- 
tude into the adiabatic regime. The present calculations, 
down to ujq = 0.01, require N gcn < 23, and could all 
be performed on small desktop computers. The valid- 
ity of our implementation was checked against calcula- 
tions from Diagrammatic Quantum Monte Carlo (data 
for the 3D polaron mass at ujq = 0.5 were provided by 
A. S. Mishchenko). 
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FIG. 1: (Color online) Polaron mass as a function of cou- 
pling for different phonon frequency, in ID (left panel) and 
3D (right panel), for cjq as indicated. 



TABLE I: Variation of numerical groundstate energy E for 



CJ 



0.1, A = 1 in ID, with the number of iterations N s . 



in the VED state selection, and corresponding number iV of 
phonon states per electron site in the Hilbert space. With 
Fock states, convergence is not achieved. With SOS, 6 digits 
of the energy are converged already for N = 10° states. The 
groundstate energy on the infinite lattice is obtained with a 
relative error of less than 10 -8 . 
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Coherent states have been used in variational^— and 
numerical approaches to the Holstcin mode l 20,42 . The 
latter two references describe a construction similar to 
ours. It appears that our construction is simpler, in that 
it avoids variational optimization of coherent state pa- 
rameters, while the additional virtues of the VED con- 
struction allow us to perform calculations for ID and 3D 
systems in the far adiabatic regime. 

III. POLARON MASS 

The fundamental quantity characterizing polaron for- 
mation is the polaron mass to defined by 



d 2 



d\k\ 2 



E(k) 



k=0 



(2) 



where E(k), the polaron dispersion, is the groundstate 
energy with momentum k. For a free electron with 
£ p = 0, the mass is mo = l/(2t). With increasing cou- 
pling, the polaron mass m grows, and the mass renor- 
malization M = m/nig becomes larger than unity. Our 
numerical method works on the infinite lattice with con- 
tinuous momentum k, which allows for a direct evaluation 
of the fc-derivative in Eq. ([2]). 

In the antiadiabatic strong coupling limit (ojq — > oo 
with g 2 fixed) the polaron mass is given by the asymp- 
totic Lang-Firsov formula M ~ cxp(g 2 ), independent of 
dimension. In contrast, the behavior in the adiabatic 
regime (ujq <C t) depends on dimension, since it is gov- 
erned by the static limit luq = 0. In the static limit the 
polaron in ID has infinite mass for all A > 0, while a 
phase transition - the self-trapping transition - from the 
free electron (M = 1) to a polaron (M = oo) occurs at a 
finite coupling strength A c rj 0.9 in 3D. Since this differ- 
ence is the most significant aspect of polaron formation 



FIG. 2: (Color online) Polaron mass as a function of phonon 
frequency for different coupling, in ID (left panel) and 3D 
(right panel), for A as indicated. 



in different dimensions, we leave out the 2D case here 
and concentrate on the ID and 3D situation. 

In Fig. [1] we show how the polaron mass at finite ujq 
interpolates between the two limits. Both for ID and 3D 
the polaron mass shows the asymptotic Lang-Firsov be- 
havior for ^> t. In ID, the curves get steeper for 
smaller ujq, while keeping their functional form. The 
curves for ujq <C t and ujq t differ only quantitatively. 
In particular, a singular value of A (or ljq) associated with 
a qualitative change does not exist. In this sense the ID 
behavior is scale-free. 

Quite different in 3D: The curves develop a sharp bend 
close to A = 1 if u) < 0.5, whose position converges to A c 
for even smaller ljq. In this case a small change of A, 
in the vicinity of A c , results in a dramatic change of to. 
Although a true phase transition does not occur— for any 
ujq > 0, already at ujq = 0.5 the behavior of the polaron 
mass is, for all practical purposes, indistinguishable from 
that for a phase transition. In the adiabatic 3D situation 
one deals either with a quasi-free electron or a very heavy 
polaron. 

Plotting the polaron mass as a function of phonon fre- 
quency (Fig. [5]) demonstrates the behavior discussed here 
in a very clear way. In ID the mass diverges for all A > 
if wo goes to zero. In 3D the mass diverges only if A > A c , 
while it approaches mo if A < A c . Just as before, now 
a small change in ujq results in a large change of m for 
A close to A c . In both situations the mass at finite fre- 
quency converges correctly to that in the static limit. 
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FIG. 3: (Color online) Polaron radius as a function of coupling 
in the static limit, in ID and 3D. Note the different scaling 
of the _R-axis on the left and right. The vertical dashed line 
indicates the critical coupling A c for the 3D adiabatic self- 
trapping transition. 
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FIG. 4: (Color online) Polaron radius as a function of cou- 
pling. Left panel: For different phonon frequencies in the 
adiabatic regime in ID. Right panel: In 3D for ljq = 0.5. The 
red dashed curves give the corresponding polaron radius in 
the static limit (from Fig. [3]), the dashed vertical line gives A c 
for the static limit. 



IV. POLARON RADIUS 

The different behavior of polaron formation in ID and 
3D, characterized through the polaron mass in the previ- 
ous section, can also be studied for the size of a polaron, 
given by the extension of the phonon cloud surrounding 
the electron. Using the electron-phonon correlation func- 
tion x(r) = SnW'oK&l+r + &n+r ) c£c„ | V>o), where IV'o) is 
the groundstate at momentum k = 0, the polaron radius 
is defined as 



R = 



E r H 2 x(r) 



1/2 



(3) 



Note that for the Holstein model, x( r ) = V e p/ W[| - 

For the static limit, the polaron radius is shown in 
Fig. [3] We observe that the radius changes continuously 
in ID. In 3D, the radius is finite only above the self- 
trapping transition at A c « 0.9, but the self-trapped po- 
laron is a very small particle, with R < 0.2. Note that in 
the static limit only the polaron radius provides informa- 
tion about the polaron properties, as the mass is infinite 
above the self-trapping transition and does not change. 

Turning to finite phonon frequencies, one should first 
understand that the polaron radius remains finite in the 
weak coupling limit. The reason is that any additional 
phonon requires a finite energy u>q, which suppresses 
phonon excitations far away from the electron where the 
energy increase is not compensated by electron-phonon 
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FIG. 5: (Color online) Electron-phonon correlation function 
x( r ) in ID for ujq — 0.01,0.04, in comparison to the static 
limit (red curve). The gray dashed curve in the left panel 
shows % f° r u o — 0.01, at weaker coupling A = 0.05. 



interaction. The limiting value of i?, as A approaches 
zero, is given by22 




(4) 



This expression implies that a large polaron can only 
exist for sufficiently small loq, e.g. R > 5 requires that 
ljq < 0.04 independent of coupling strength. A different 
question is the existence of a large and heavy polaron. As 
the results for the static limit indicate, such a polaron can 
only exist in ID, at very small ujq <C t and small coupling 
A < 0.5. 

In Fig.0]we show the polaron radius in ID and 3D. One 
should note that the evaluation of Eq. is not possible 
with extremely high accuracy, due to the summation in 
the numerator containing |r| 2 . In fact, the expression 
is not even well conditioned since states with arbitrarily 
small energy difference can have totally different R (just 
excite a phonon very far from the electron). For the 
present purpose, we calculate R under the assumption 
that x( r ) decays exponentially for sufficiently large |r|, 
and the resulting error is about 3% in the most difficult 
situation (u>o = 0.01, A ~ 0.5). 

As expected, R <C 1 in 3D if A > A c . In ID we find 
that, at the smallest frequency ujq = 0.01 studied, a po- 
laron with radius R > 2 (R > 3) exists for A < 0.43 
(A < 0.3). If we compare these values with the polaron 
mass in Fig. [TJ we find, e.g. for A = 0.25, a large and 
heavy polaron with R > 4, M > 10. 

To support our claim that polarons for these parame- 
ters correspond to what is often called an adiabatic po- 
laron in the literature we show in Fig. [5] the electron- 
phonon correlation function x( r ) m ID for loq = 0.01, 
A = 0.25 in comparison to that for the static limit, which 
closely resemble each other. The agreement becomes 
even better for A = 0.4, although then the polaron with 
R > 2 is not really large. Note that in order to find such 
agreement one has to use phonon frequencies as small as 
ujq = 0.01. Already at u>q = 0.04 the curves do not fit. 
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FIG. 6: (Color online) Curves of constant mass, and curves of 
constant radius in ID. The left panel is a magnification of the 
adiabatic regime cjo < 0.25, which corresponds to the small 
orange rectangle in the right panel. 



V. SURVEY OF THE ID POLARON 




0.25 2 4 6 8 10 

FIG. 7: (Color online) Overview of different polaron regimes 
in ID. Only for the notation in this figure, we call a polaron 
light (heavy) if M < 2 (M > 10), and large (small) if R > 3 
(7? < 0.5). With this terminology, heavy large polarons are 
found in the region "HL". The gray dashed rectangle gives 
the parameter region shown in the left panel in Fig. [6] For 
better visibility, the region of smaller ujo , A values is magnified 
through a continuous variation of the scaling along each axis. 



VI. OPTICAL RESPONSE IN ID 



We have seen in the previous two sections that in 
the adiabatic regime in 3D light polarons evolve into 
very heavy and very small polarons through a rapid 
crossover. Only in ID a continuous renormalization of 
the polaron mass and radius occurs in an extended re- 
gion of parameter space. A map of that region is given 
in Fig. [S] where we plot curves of constant mass and of 
constant radius. If we assign names to certain regimes 
bordered by the M = 2, 10 and R = 0.5,3 curves, we 
get Fig. which indicates the tendencies in the evolu- 
tion from a light to a heavy, or a large to a small po- 
laron with changing model parameters. Independent of 
the precise choice of wording, large and heavy polarons 
exist for sufficiently small loq. The existence of such po- 
larons is specific for ID, and a similar regime does not 
appear in 3D, where heavy polarons are always small as 
we discussed in the previous section. Even in ID, it re- 
quires very small phonon frequencies. In combination 
with the comparison of the electron-phonon correlation 
function in Fig. [5] we can therefore refine the conven- 
tional adiabaticity condition ujq < t to a much stricter 
statement, namely wo/t <C 0.1 or uio/t < 0.01. The spe- 
cific ID physics of large heavy polarons thus occurs only 
in a tiny region of parameter space. Turning to larger 
phonon frequencies, the polaron is always small. For 
u>o/t 3> 1, the polaron is almost point-like. Rearrang- 
ing Eq. (J3J) gives x(0)/S n x( n ) > 1 - R 2 , so that e.g. 
for R = 0.5 at least 75% of the total lattice displace- 
ment are located at the electron. In this limit, the only 
significant polaron parameter is the average number of 
phonons g 2 = e p /uio - instead of the two parameters ojq 
and A -, the mass is given by the asymptotic Lang-Firsov 
expression M — exp(g 2 ), and the differences between di- 
mensions vanish. 



The physics of adiabatic, or large and heavy, polarons 
is essentially different from that of small polarons in the 
antiadiabatic strong coupling regime. The difference ex- 
tends beyond the differences in mass or size of the po- 
laron discussed so far, but affects its very structure, and 
thus relates to the question in which sense a polaron can 
be understood as an itinerant quasiparticl o 22 ' 44 ' 45 . For 
a characterization of the quasiparticle properties we can 
study its response to electric fields. The (regular part of 
the) optical conductivity o- reg (u;) is defined through 

a rcg H = 7T V ( ^' Wo)|2 [5(w - ui) + 5(uj + cut)} , (5) 

— ' UJ] 

l>0 

with the current operator j = it J2 n c n+i c n ~ c n c n+i an d 
the eigenstates with excitation energy w/. Further 
information is contained in the f-sum rule 

2irV + 2 / a reg (u)du = E kin , (6) 
Jo 

where E kin = -t(tp \ E n c n+i c n + clcn+i^o) is the ki- 
netic energy. As a basic transport quantity here the 
Drudc weight T> appears, which is related to the polaron 
mass m by T> = l/(2m). 

The reader should be aware that in a numerical cal- 
culation of (T rC g(w) it is impossible to capture the exact 
position and weight of all excited states. For our results, 
the resulting error in integrated quantities, e.g. the first 
few moments of <7 reg (u;), is negligible. For cr rGg (u;) itself 
a compromise between finite spectral resolution and the 
finite level spacing in the numerical Hilbert space has to 
be found. Consequently, the shape, envelope and peak 
positions in the pictures for a Teg (uj) (see Fig. |H] below) 



FIG. 8: (Color online) Depiction of the f-sum rule, with 
kinetic energy _Ekin (solid), integrated regular conductivity 
SVcg = J °° o" reg (u/)du/ (dashed), and Drude weight T> (dash- 
dotted) in ID. The long-dashed red curve in the left column 
gives the kinetic energy in the static limit for comparison. For 
ujo = 0.01, it lies on top of the .Ekm-curve. 



are correct, but it should be clear that tiny spikes on top 
of the curves are a kind of numerical noise. 

We show the evaluation of the f-sum rule in Fig. [8] 
Note that, working on an infinite lattice, the f-sum rule 
is satisfied in all parameter regimes. In particular, the 
Drude weight obtained either from the f-sum rule via 
cr reg (w) or from the polaron mass differs by less than 
10" 4 . 

As we have discussed above, the polaron in the an- 
tiadiabatic strong coupling regime with R < 1 attains 
a large mass renormalization only through local phonon 
excitations tightly bound to the electron. Such a polaron 
is always a very small particle, with strong correlations 
between electron and phonon dynamics. We therefore 
expect that, with increasing coupling, the kinetic en- 
ergy and Drude weight simultaneously become small as 
the polaron gets less mobile, and that at the same time 
the optical conductivity decreases since a very small ob- 
ject is less susceptible to optical fields. As this expecta- 
tion is confirmed in Fig. [8] (right bottom panel), we can 
say that in this type of polaron the electron is dressed 
with phonons, resulting in a new polaronic quasiparticle 
emerging as a joint electron-phonon entity. 

The physics is essentially different in the adiabatic 
regime of heavy but large polarons. For such adiabatic 
polarons electron and phonon dynamics partially decou- 
ple, anticipating the static limit where the electron moves 
in the rigid phonon configuration of the lattice displace- 
ment. That there translational symmetry is broken and 
the polaron has infinite mass is equivalent to the fact that 
states of a static oscillator with different position, and 
vanishing width, have zero overlap. At small but finite 
wo the translated phonon configurations have small over- 
lap, and the polaron acquires a finite but eventually large 
mass even for large radius. That the radius remains large 
means that in contrast to the antiadiabatic situation the 



FIG. 9: (Color online) Optical conductivity a rcg (uj) in ID. 
Left column: For small u>o. The red dashed curves in the 
upper two panels give the conductivity in the static limit 
ujo = 0, with same A. Right column: For large u>o- The 
red dashed curves show the approximate Gaussian expression 
for the conductivity obtained in small-polaron theor y 46 ' 47 . 
The insets compare the integrated conductivity S rcs (co) = 
J o~ IC g(u>')du)' (in arbitrary units). 



mass increase is not so much a result of dressing the elec- 
tron with phonons, but rather of the immobility of the 
extended lattice displacement as a whole. The electron 
motion relative to the lattice displacement is less sup- 
pressed, which results in the large radius. Again, we see 
this behavior in the f-sum rule in Fig. |8] (left top panel). 
Now, even if the Drude weight close to zero indicates the 
immobility of the polaron, the kinetic energy and opti- 
cal conductivity remain large. The opposite behavior of 
these quantities can be attributed to the different con- 
stituents of the polaron. The Drude weight is small since 
the phonon part suppresses motion of the whole polaron, 
but the kinetic energy is large due to the electron motion 
relative to the phonons. Transitions between different 
electronic states in the almost rigid phonon background 
give a large contribution to the optical conductivity. If 
electron and phonon motion decouple in this way an adi- 
abatic polaron can be understood as a quasiparticle only 
with certain restrictions. The polaron is rather a compos- 
ite electron-phonon object, where electron and phonons 
influence each other but remain separable, instead of that 
joint electron-phonon entity found in the antiadiabatic 
regime where the electron and the phonons lose their in- 
dependence in favor of a new quasiparticle. As the f-sum 
rule shows this difference persists through the entire cou- 
pling range. 

Differences between the adiabatic and antiadiabatic 
regime appear also in the optical conductivity (see 
Fig. [5]). At large phonon frequencies and strong coupling 
the optical conductivity is dominated by phonon-assisted 
transitions to excited polaron states, which change the 
average phonon number by an integer. The single peaks 
are thus separated by ujq, and resolve the Poisson dis- 
tribution of phonons in the polaron groundstate. The 
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skewed Gaussian envelope of a Teg (uj) is a signature of 
antiadiabatic small polaron a 46 ' 47 . Note that the inte- 
grated optical conductivity <SV eg is small in this regime. 
For intermediate phonon frequencies complicated struc- 
tures develop, with a contribution from different polaron 
bands whose width is no longer small compared to ljq. 
In the far adiabatic regime the basic mechanism for op- 
tical absorption has changed. It is now dominated by 
electronic transitions in the potential well generated by 
the lattice displacement. Of course, a phonon admixture 
is present for any finite lu . Only in the static limit the 
phonon configuration remains unchanged during an elec- 
tronic transition, and the Franck-Condon effect is fully 
realized. 

To summarize, the change in the optical conductiv- 
ity with phonon frequency reflects very clearly the two 
different time scales realized in the polaron problem. De- 
pending on whether the electron or the phonons set the 
fastest time scale, we observe different physics. In the an- 
tiadiabatic regime we find a well-defined polaronic quasi- 
particle whose behavior is governed by local phonon ex- 
citations dressing the electron and the arising strong cor- 
relations between electron and phonon dynamics. Or we 
find, in the adiabatic regime, an electron moving in a 
phononic background whose properties depend only on 
the average electron motion, while stronger correlations 
are missing. These differences and their consequences for 
polaron transport and dynamics have been studied here 
in terms of the optical conductivity and f-sum rule. 



VII. CONCLUSION AND OUTLOOK 

A description of lattice polaron formation requires the 
treatment of the two different time-scales of electron and 
phonon dynamics, their mutual influence and accompa- 
nying renormalization, and the non-trivial quantum cor- 
relations that develop. The adiabatic regime is especially 
difficult, which is typical for many problems where one 
tries to perform the classical limit for some part of a 
system while keeping the full quantum dynamics of the 
remainder. 

For the Holstein model one is in the fortunate posi- 
tion that physical concepts or approximate theories can 
be checked against very accurate data provided by ad- 
vanced numerical calculations. In that way the valid- 
ity of important ideas, potentially with generalizations 
to less tractable situations, can be assessed, and some 
of the controversies in the polaron literature may be re- 
solved. In the present contribution we develop our phys- 
ical discussion along the lines indicated by such accurate 
data. The new quality of this contribution is our abil- 
ity to cover the entire coupling and frequency range, in 
particular the far adiabatic regime including heavy large 
polarons in ID. This progress becomes possible by the 
new numerical method introduced here. 

Our results support a particular picture of polaron for- 
mation, which is in agreement only with parts of the liter- 



ature. Let us summarize some important findings. (I) Al- 
though all groundstate quantities vary continuously with 
coupling and phonon frequency, the steady change in ID 
is replaced with a rapid crossover in the adiabatic regime 
in 3D, which closely resembles a true phase transition. 
While no characteristic energy scale exists in ID, it is 
set in 3D by the critical coupling A c of the self-trapping 
transition in the static limit. 

(II) All quantities converge (or diverge) to the static limit 
in a simple and clear way. Nevertheless quantitative 
agreement is obtained only for very small phonon fre- 
quencies OJa/t <C 0.1. 

(III) Only for such small phonon frequencies heavy and 
large polarons occur in ID (e.g. luq < 0.02 for M > 5 
and R > 4). A large radius requires a small phonon fre- 
quency in any case. 

(IV) Small and large polarons have significantly different 
transport and optical properties. A large and heavy po- 
laron is characterized through a small Drude weight but 
large optical absorption. A small and heavy polaron has 
small Drude weight and small optical absorption. 

(V) The small-polaron approximation for the optical con- 
ductivity is strictly valid only for very large phonon fre- 
quencies, while in general pronounced structures distort 
the Gaussian envelope. The overall shape of the envelope 
however remains comparable also for smaller phonon fre- 
quencies. 

These findings suggest certain physical interpretations. 
First, with all necessary caution, it may make sense to 
think of a self-trapping transition in 3D even at finite 
phonon frequency. There is nothing comparable in ID. 
Second, the notion of an itinerant polaronic quasiparticle 
is to be used only with care: The nature of the quasi- 
particle differs at large and small phonon frequency (cf. 
Sec. IVip . Third, Holstein polarons tend to be small ob- 
jects. This should largely prevent polaron dissociation, 
which may occur in a kind of polaronic Mott transition 
if the phonon clouds overlap at finite polaron density. It 
also prevents instabilities in anisotropic systems with a 
preferred direction of motion, where only the large adia- 
batic quasi-lD polarons arc destabilized by electron hop- 
ping perpendicular to this directio n 48 ' 49 . 

In principle, our results allow for direct comparison 
with, and hence validation of, approximate polaron the- 
ories. We have not tried to do so here, since we want 
our results to tell their own tale. Our discussion should 
have shown that the numerical data itself do not prevent 
a physical interpretation. But the interpretation adheres 
to direct terms - the polaron is heavy or light, small or 
large -, and simple pictures - such as a phonon-dressed 
electron vs. an electron moving in the potential gener- 
ated by the lattice displacement. There is no necessity to 
invoke complicated constructions which do not find their 
counterpart in the few quantities truly characterizing the 
polaron. 

The present problem, i.e. the single Holstein polaron at 
zero temperature, is treated here exhaustively Of course, 
it is possible to go to even smaller phonon frequencies, 
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larger couplings or other more extreme situations, but 
new physical behavior will not occur. There are how- 
ever many other topics in polaron physics that remain 
unsettled. With a few singular exceptions none of them 
have been accessed with the same rigor as the present 
problem. Within the framework of the Holstein model, 
transport at finite temperature and collective effects at fi- 
nite density are probably the two most important issues. 
Proceeding beyond the Holstein model, in particular on 
trying to extend the discussion to polarons in novel ma- 
terials such as graphene 5 ^, new questions emerge. These 
include polaron formation in the presence of impurities^ 
and disorde r 51 ' 52 , the influence of phonon dispersion, of 
different types of electron-phonon coupling, of lattice ge- 
ometry anisotropy^ or confinement, or the dynamics of 
polaron formation^ out of equilibrium, realized e.g. in 
quantum dots or molecular aggregates. Some of these 
questions are directly connected with the present study. 
Does, e.g., longer-ranged electron-phonon interaction fa- 
vor the formation of large and heavy polarons? Although 
the polaron becomes larger then, its mass will decrease, 
and the existence of a large heavy polaron with corre- 
sponding optical signatures may crucially depend on a 
number of material parameters. The classification of 
these parameters and of their influence on polaron prop- 
erties is central to giving studies such as the present one 
relevance beyond the rigorous treatment of a model. To 
make such investigations possible is certainly the moti- 
vation to develop powerful high-precision numerics for 
the polaron problem, while their application provides us 
with the necessary information to understand real po- 
larons and the models to describe them. 
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Appendix A: Shifted Oscillator States 

The SOS \n) a , with n > 0, form a Hilbert space basis 
for a single harmonic oscillator. For a ^ 0, they are 
characterized through two properties: (i) orthonormality 
(n\m) a = 5 nm , (ii) the n-th SOS is a linear combination 



n 
m=0 



(Al) 



These conditions imply that the SOS can be con- 
structed through Gram-Schmidt-orthogonalization of the 
equally spaced coherent states \ma) c , in the order of in- 
creasing m = 0,1,2,.. .,oo. Moreover, the operator b is 
given by a upper triangular matrix in the SOS basis. It is 
possible to obtain explicit expressions for the coefficients 



,(n) 



na- 

fc=i 



2A'- 



1/2 



(A2) 



where a = e and the q-binomials [m] q (f° r a defini- 
tion, see e.g. Ref. l54h are polynomials in q. All relevant 
matrix elements can also be computed with a numerical 
Gram-Schmidt procedure. 

The above definition applies to all a ^ 0. The limit 
d — y exists, and the SOS reproduce the Fock states in 
this limit. In the opposite limit a — > oo, they become 
asymptotically equal to the coherent states \ma) c , since 
the overlap between different coherent states is small. 

From the explicit expressions for c"m we find that 
(jilfol?!)^ = ?ier. For real er, it follows that the oscilla- 
tor elongation in the n-th SOS is given by 2na. Thus, 
SOS for a > (cr < 0) should be used for states with 
positive (negative) elongation. 

To quantify the benefit of SOS over Fock states for 
groundstate calculations, we consider a shifted oscillator 
with Hamiltonian 



H osc = to tfb - gtd (tf + b) 



(A3) 



of the n + 1 coherent states \ma) Cl < m < n. In 
particular, |0) CT = |0) c = |vac). 



The groundstate of this Hamiltonian, with energy — e p = 
— g 2 uiQ, is a coherent state \g) Cl which contains g 2 bosons 
on average. In a numerical calculation an approximate 
groundstate of H osc is calculated in a truncated Hilbert 
space, e.g. consisting of the first N + 1 Fock states. The 
corresponding error in the groundstate energy is shown 
in Fig. [TDI As expected, the error is large when TV < g 2 . 

Using SOS instead of Fock states, the error depends on 
cr. In the present example, the error is zero for a = g/N , 
since the coherent state \g) c is & linear combination of the 
first two SOS. For a meaningful error estimate we there- 
fore fix a — g ma,x /N for the maximal g max considered. 
The virtue of the SOS is that the error is small also for 
all intermediate < g < ,g max (see Fig. [TO]). Note that 
the error drops to zero at those values of g for which the 
true groundstate is contained in the truncated Hilbert 
space spanned by the N +1 SOS. 

We also see how the SOS prefer positive elongations 
for cr > 0, and the error becomes large if a and g have 
opposite sign. In a real calculation, the variational de- 
termination of the parameter a would prevent this situ- 
ation. For the Holstein model ([I}, only positive elonga- 
tions occurs and always a > 0. Note that in the worst 
case the optimal a is zero, and the Fock states are re- 
covered. Therefore, the introduction of SOS can only 
improve groundstate results obtained with Fock states. 
In the given example, the error with 10 SOS is one order 
of magnitude smaller than that with 100 Fock states. 
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FIG. 10: Absolute error ol the groundstate energy of the 
Hamiltonian ()A3[) with ujq = 0.01, for a maximal g 2 = 100, 
as explained in the text. Left panel: Error using N + 1 
Fock states. Right panel: Error using N = 10 SOS with 
a = g max /N = 1, as a function of g. 

The SOS construction can be modified in various ways. 
For example, if oscillator states with large positive and 
negative elongations should be represented, one can ap- 
ply a global shift S(£), constructing the SOS from coher- 
ent states |m<7+£) c . The additional parameter £ provides 
the necessary freedom to deal with such a situation. 

We introduced the SOS mainly for oscillator ground- 
states. They also work for excited states, but we should 
distinguish two cases. First, for a single excited state, the 
SOS with variational determination of a will succeed and 



can provide us with good results even in cases when Fock 
states fail. Or, second, we may need all excited states at 
once, such as in the calculation of the optical conductivity 
Creg(uj). In this case, the fundamental problem of any ED 
calculation is that a finite Hilbert space does not allow 
for the representation of infinitely many excited states. 
No basis construction can overcome this obstruction. It 
is however possible to include sufficiently many states in 
the Hilbert space, such that the Hilbert space truncation 
shows up only through a kind of numerical noise, while 
the relevant physical structures are resolved. For exam- 
ple, the equally spaced peaks in the lower right panel in 
Fig. [9] are physically meaningful, while the tiny spikes 
in the middle right panel are numerical noise. The SOS 
construction guarantees the correct envelope and shape 
of cr r eg(w). For example, in the middle left panel in Fig. [HI 
a calculation with Fock states instead of SOS results in 
a curve that is shifted to smaller energies. Again, the 
Fock states cannot account for the substantial oscillator 
elongation in the groundstate and excited states in the 
adiabatic regime. To summarize, integrated or averaged 
properties of excited states are obtained correctly with 
SOS, but not with Fock states. For correct spectral prop- 
erties, e.g. peak positions at high resolution, additional 
effort is needed in any case. 
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